## How do people evaluate foreign aid to “nasty” regimes?
## Tobias Heinrich & Yoshiharu Kobayashi
## British Journal of Political Science
#########################################################


## Run regressions on expertimental treatments


## Load data
############
load("output/Modified MTurk Data.Rdata")
load("output/Weight vectors.Rdata")


## Estimate Model for experimental variables
############################################
mod_form <- list(A1=Rating ~ 1 + I(C == 25) + I(C == 75) + I(B == 1) + I(B == 2) + I(B == 3) + I(B == 4) +
                   + I(PI == 1) + I(PI == 2) + I(PI == 3) + I(PI == 4) + I(PI == 5) +
                   + I(PI == 2 & RA == 0) + I(PI == 3 & RA == 0) + I(PI == 4 & RA == 0) + I(PI == 5 & RA == 0)  +
                   + I(I(PI == 2 & RT == 1) * RA) + I(I(PI == 2 & RT == 2) * RA) + I(I(PI == 2 & RT == 3) * RA) +
                   + I(I(PI == 3 & RT == 1) * RA) + I(I(PI == 3 & RT == 2) * RA) + I(I(PI == 3 & RT == 3) * RA) +
                   + I(I(PI == 4 & RT == 1) * RA) + I(I(PI == 4 & RT == 2) * RA) + I(I(PI == 4 & RT == 3) * RA) +
                   + I(I(PI == 5 & RT == 1) * RA) + I(I(PI == 5 & RT == 2) * RA) + I(I(PI == 5 & RT == 3) * RA),
                 A2=Rating ~ 1 + I(C == 25) + I(C == 75) + I(B == 1) + I(B == 2) + I(B == 3) + I(B == 4) + I(PI == 1) + I(PI == 2) + I(PI == 3) + 
                   I(PI == 4) + I(PI == 5),
                 A3=Rating ~ 1 + I(C == 25) + I(C == 75) + I(B == 1)  + I(B == 2)  + I(B == 3)  + I(B == 4) +
                   + I(PI == 1) + I(PI == 2) + I(PI == 3) + I(PI == 4) + I(PI == 5) +
                   + I(PI == 2 & B == 1) + I(PI == 3 & B == 1) + I(PI == 4 & B == 1) + I(PI == 5 & B == 1) +
                   + I(PI == 2 & B == 2) + I(PI == 3 & B == 2) + I(PI == 4 & B == 2) + I(PI == 5 & B == 2) +
                   + I(PI == 2 & B == 3) + I(PI == 3 & B == 3) + I(PI == 4 & B == 3) + I(PI == 5 & B == 3) +
                   + I(PI == 2 & B == 4)  + I(PI == 3 & B == 4) + I(PI == 4 & B == 4) + I(PI == 5 & B == 4),
                 A4=Rating ~ 1 + I(C == 25) + I(C == 75) + I(B == 1) + I(B == 2) + I(B == 3) + I(B == 4) +
                   + I(PI == 1) + I(PI == 2) + I(PI == 3) + I(PI == 4) + I(PI == 5) +  
                   + I(PI == 2 & C == 25) + I(PI == 3 & C == 25) + I(PI == 4 & C == 25) + I(PI == 5 & C == 25) +
                   + I(PI == 2 & C == 75) + I(PI == 3 & C == 75) + I(PI == 4 & C == 75) + I(PI == 5 & C == 75))


all_exp_mod <- vector("list", length(mod_form))

## Run the regressions
for(i in 1:length(mod_form))
{
  all_exp_mod[[i]] <- vector("list", ncol(weight_vec[[1]]))
  for(j in 1:ncol(weight_vec[[1]]))
  {
    if(names(mod_form)[i] == "A1")
    {
      dat <- data
      w <- weight_vec[[1]][,j]
    }
    if(names(mod_form)[i] != "A1")
    {
      dat <- subset(data, RA == 0)
      w <- weight_vec[[2]][,j]
    }
    
    
    ## Estimate model and calculate cluster-BS adjusted VCO 
    mod <- lm(formula=as.formula(mod_form[[i]]), data=dat, weight=w)              
    cluster <- dat$ID
    clusters <- unique(cluster)
    nc <- length(clusters)
    Obsno <- split(1:nrow(dat), cluster)
    nvars <- length(coef(mod))
    bar <- rep(0, nvars)
    cov <- matrix(0, nvars, nvars)    
    for(iter in 1:n_bs)
    {
      ## Draw clusters
      k <- sample(clusters, nc, replace = TRUE)
      obs <- unlist(Obsno[k])
      
      mod0 <- lm(formula=as.formula(mod_form[[i]]), data=dat[obs,], weights=w[obs])
      bar <- bar + coef(mod0)
      cov <- cov + coef(mod0) %*% t(coef(mod0))
    }    
    ## Obtain cluster-adjusted VCOV matrix
    bar <- bar/n_bs
    cov <- (cov - n_bs * bar %*% t(bar))/(n_bs - 1)  
    
    all_exp_mod[[i]][[j]]$Model <- mod
    all_exp_mod[[i]][[j]]$BS_VCV <- cov
    all_exp_mod[[i]][[j]]$TheForm <- as.formula(mod_form[[i]])
    
    
    ## Calculate quantities for the effects graphs
    ##############################################
    
    ####################
    ## Analysis for "A1"
    ## PI x RA
    ####################
    if(names(mod_form)[i] == "A1")
    {
      ## Conditional effects perspective
      ###############################
      beta <- rmvnorm(1000, coef(mod), cov)
      tmp_placebo <- beta[, "I(PI == 1)TRUE"]
      output <- data.frame(Case="Placebo", Which="None", Info="A1", 
                           q_l=quantile(beta[, "I(PI == 1)TRUE"], .025),
                           q_u=quantile(beta[, "I(PI == 1)TRUE"], .975),
                           q_m=quantile(beta[, "I(PI == 1)TRUE"], .5),
                           m=mean(beta[, "I(PI == 1)TRUE"]),
                           sd=sd(beta[, "I(PI == 1)TRUE"]),
                           d_means=NA, d_means_ci_l=NA, d_means_ci_h=NA,
                           ratio_placebo=1, ratio_placebo_ci_l=1,
                           ratio_placebo_ci_h=1)
      
      ra_levs <- seq(from=0, to=25, by=1)
      for(k in 1:4)
      {
        ## Calculate no response case
        tmp_base <- rowSums(beta[, c(paste0("I(PI == ", k+1, ")TRUE"),
                                     paste0("I(PI == ", k+1, " & RA == 0)TRUE"))])
        output <- rbind(data.frame(Case=k + 1, Which="No remedial aid", Info="A1", q_l=quantile(tmp_base, .025),
                                   q_u=quantile(tmp_base , .975), q_m=quantile(tmp_base, .5),
                                   m=mean(tmp_base), sd=sd(tmp_base), d_means=NA, d_means_ci_l=NA, d_means_ci_h=NA,
                                   ratio_placebo=mean(tmp_base/tmp_placebo),
                                   ratio_placebo_ci_l=quantile(tmp_base/tmp_placebo, .025),
                                   ratio_placebo_ci_h=quantile(tmp_base/tmp_placebo, .975)),
                        output)
        
        ## Calculate best response's conditional effects
        out <- array(NA, dim=c(4, length(ra_levs), 1000))
        for(m in 1:4)
        {
          for(n in 1:length(ra_levs))
          {
            out[m, n, ] <- beta[, paste0("I(PI == ", k+1, ")TRUE")] 
            if(m > 1) out[m, n, ] <- out[m, n, ] + beta[, paste0("I(I(PI == ", k+1, " & RT == ", m-1, ") * RA)")] * ra_levs[n]
          }
        }
        
        tmp <- aaply(.data=out, .margins=3, .fun=function(x) max(c(x)))
        output <- rbind(data.frame(Case=k + 1, Which="Best response", Info="A1", q_l=quantile(tmp, .025),
                                   q_u=quantile(tmp, .975), q_m=median(tmp), m=mean(tmp), sd=sd(tmp),
                                   d_means=mean(tmp)-mean(tmp_base), d_means_ci_l=quantile(tmp-tmp_base, .025),
                                   d_means_ci_h=quantile(tmp-tmp_base, .975),
                                   ratio_placebo=mean(tmp/tmp_placebo),
                                   ratio_placebo_ci_l=quantile(tmp/tmp_placebo, .025),
                                   ratio_placebo_ci_h=quantile(tmp/tmp_placebo, .975)),
                        output)
      }
    }
    
    
    ####################
    ## Analysis for "A2"
    ## Unconditional effects
    ####################
    if(names(mod_form)[i] == "A2")
    {
      output <- c()
    }
    
    ####################
    ## Analysis for "A3"
    ## PI x B
    ####################
    if(names(mod_form)[i] == "A3")
    {
      ## Conditional effects perspective
      ##################################
      beta <- rmvnorm(1000, coef(mod), cov)
      tmp_placebo <- beta[, "I(PI == 1)TRUE"]
      output <- data.frame(Case="Placebo", Which="Baseline", Info="A3", 
                           q_l=quantile(beta[, "I(PI == 1)TRUE"], .025),
                           q_u=quantile(beta[, "I(PI == 1)TRUE"], .975),
                           q_m=quantile(beta[, "I(PI == 1)TRUE"], .5),
                           m=mean(beta[, "I(PI == 1)TRUE"]),
                           sd=sd(beta[, "I(PI == 1)TRUE"]),
                           d_means=NA, d_means_ci_l=NA, d_means_ci_h=NA,
                           ratio_placebo=1, ratio_placebo_ci_l=1,
                           ratio_placebo_ci_h=1)
      
      for(k in 1:4)
      {
        ## Calculate effect PI | B=0
        tmp_base <- beta[, paste0("I(PI == ", k+1, ")TRUE")]
        output <- rbind(data.frame(Case=k + 1, Which="Baseline", Info="A3", q_l=quantile(tmp_base, .025),
                                   q_u=quantile(tmp_base, .975), q_m=quantile(tmp_base, .5), m=mean(tmp_base), sd=sd(tmp_base),
                                   d_means=NA, d_means_ci_l=NA, d_means_ci_h=NA,
                                   ratio_placebo=mean(tmp_base/tmp_placebo),
                                   ratio_placebo_ci_l=quantile(tmp_base/tmp_placebo, .025),
                                   ratio_placebo_ci_h=quantile(tmp_base/tmp_placebo, .975)),
                        output)
        
        ## Calculate effect PI | B = m
        for(m in 1:4)
        {
          tmp <- rowSums(beta[, c(paste0("I(PI == ", k+1, ")TRUE"),
                                  paste0("I(PI == ", k+1, " & B == ", m, ")TRUE"))])
          output <- rbind(data.frame(Case=k + 1, Which=paste0("Benefits", m), Info="A3", q_l=quantile(tmp, .025),
                                     q_u=quantile(tmp, .975), q_m=quantile(tmp, .5), m=mean(tmp), sd=sd(tmp),
                                     d_means=mean(tmp) - mean(tmp_base),
                                     d_means_ci_l=quantile(tmp - tmp_base, .025), 
                                     d_means_ci_h=quantile(tmp - tmp_base, .975),
                                     ratio_placebo=mean(tmp/tmp_placebo),
                                     ratio_placebo_ci_l=quantile(tmp/tmp_placebo, .025),
                                     ratio_placebo_ci_h=quantile(tmp/tmp_placebo, .975)),
                          output)
        }    
      }
    }
    
    ####################
    ## Analysis for "A4"
    ## PI x C
    ####################
    if(names(mod_form)[i] == "A4")
    {
      ## Conditional effects perspective
      ##################################
      beta <- rmvnorm(1000, coef(mod), cov)
      tmp_placebo <- beta[, "I(PI == 1)TRUE"]
      output <- data.frame(Case="Placebo", Which="None", Info="A4", 
                           q_l=quantile(beta[, "I(PI == 1)TRUE"], .025),
                           q_u=quantile(beta[, "I(PI == 1)TRUE"], .975),
                           q_m=quantile(beta[, "I(PI == 1)TRUE"], .5),
                           m=mean(beta[, "I(PI == 1)TRUE"]),
                           sd=sd(beta[, "I(PI == 1)TRUE"]),
                           d_means=NA, d_means_ci_l=NA, d_means_ci_h=NA,
                           ratio_placebo=1, ratio_placebo_ci_l=1,
                           ratio_placebo_ci_h=1)
      
      for(k in 1:4)
      {
        ## Calculate effect PI | C=75
        tmp_75 <- rowSums(beta[, c(paste0("I(PI == ", k+1, ")TRUE"),
                                paste0("I(PI == ", k+1, " & C == 75)TRUE"))])
        output <- rbind(data.frame(Case=k + 1, Which=75, Info="A4", q_l=quantile(tmp_75, .025), 
                                   q_u=quantile(tmp_75, .975), q_m=quantile(tmp_75, .5),
                                   m=mean(tmp_75), sd=sd(tmp_75),
                                   d_means=NA, d_means_ci_l=NA, 
                                   d_means_ci_h=NA,
                                   ratio_placebo=mean(tmp_75/tmp_placebo),
                                   ratio_placebo_ci_l=quantile(tmp_75/tmp_placebo, .025),
                                   ratio_placebo_ci_h=quantile(tmp_75/tmp_placebo, .975)),
                        output)

        ## Calculate effect PI | C=25
        tmp <- rowSums(beta[, c(paste0("I(PI == ", k+1, ")TRUE"),
                                paste0("I(PI == ", k+1, " & C == 25)TRUE"))])
        output <- rbind(data.frame(Case=k + 1, Which=25, Info="A4", q_l=quantile(tmp, .025), 
                                   q_u=quantile(tmp, .975), q_m=quantile(tmp, .5), m=mean(tmp), sd=sd(tmp),
                                   d_means=mean(tmp) - mean(tmp_75),
                                   d_means_ci_l=quantile(tmp - tmp_75, .025), 
                                   d_means_ci_h=quantile(tmp - tmp_75, .975),
                                   ratio_placebo=mean(tmp/tmp_placebo),
                                   ratio_placebo_ci_l=quantile(tmp/tmp_placebo, .025),
                                   ratio_placebo_ci_h=quantile(tmp/tmp_placebo, .975)),
                        output)
        
        ## Calculate effect PI | C=50
        tmp <- beta[, paste0("I(PI == ", k+1, ")TRUE")]
        output <- rbind(data.frame(Case=k + 1, Which=50, Info="A4", q_l=quantile(tmp, .025),
                                   q_u=quantile(tmp, .975), q_m=quantile(tmp, .5),
                                   m=mean(tmp), sd=sd(tmp),
                                   d_means=mean(tmp) - mean(tmp_75),
                                   d_means_ci_l=quantile(tmp - tmp_75, .025), 
                                   d_means_ci_h=quantile(tmp - tmp_75, .975),
                                   ratio_placebo=mean(tmp/tmp_placebo),
                                   ratio_placebo_ci_l=quantile(tmp/tmp_placebo, .025),
                                   ratio_placebo_ci_h=quantile(tmp/tmp_placebo, .975)),
                        output)
        
      }
    }
    
    output$w <- j
    output$What <- "Placebo"
    output$What[output$Case == 5] <- "Media\ncrackdown"
    output$What[output$Case == 4] <- "Torture"
    output$What[output$Case == 3] <- "Rigged\nelections"
    output$What[output$Case == 2] <- "Aid\ntheft"
    output$What <- factor(output$What, levels=c("Placebo", "Aid\ntheft", "Rigged\nelections",
                                                "Media\ncrackdown", "Torture"))

    all_exp_mod[[i]][[j]]$BestResponseMR <- output
    
    ## Format and save output to mention in the text
    round_vars <- c("q_l", "q_u", "q_m", "m", "sd", "d_means", "d_means_ci_l", "d_means_ci_h",
                    "ratio_placebo", "ratio_placebo_ci_l", "ratio_placebo_ci_h")
    if(i != 2) for(k in 1:length(round_vars)) output[, round_vars[k]] <- round(output[, round_vars[k]], di=2)    
    write.csv(x=output, file=paste0(getwd(), "/output/tables/Output-Case", i, "-W", j, ".csv"))
  }
}


## Save files
save(all_exp_mod, file=paste0(getwd(), "/output/ExpModels.Rdata"))

